% code used to create Figs. 6 in "An alternative pathway for delivery of 
% power by outer hair cells to cochlear traveling waves", Samaras and Meaud
% plot  tonic displacement figur
% Julien Meaud
% julien.meaud@me.gatech.edu
% ------------------------------------------------------------------------

close all
clear all

setPlottingOptions

normalizedTime=1;
nOffset=17.5;  %number of cycle offset 

FS=10;  % fontsize

Nmodel=2; % 2 to plot M1 and M4; 3 to plot M1, M4, and M5

for iModel=1:Nmodel
    if iModel==1    %M1   
          var=load('TimeDomain_CdcCohc/TB_G35_Act1_85B1_15A_rEps3_1_rKdcKohc_1000_Krl_BL_Cdc0_0_25_ohc0_0_25_80PT20000.mat');       
    elseif iModel==2 %M4
        var=load('TimeDomain_CdcCohc/TB_G35_Act3_0B1_15A_rEps3_1_rKdcKohc_5_Krl_90BL_Cdc0_0_25_Cohc0_0_25_80PT20000.mat');
    elseif iModel==3 %M5       
        var=load('TimeDomain_CdcCohc/TB_G35_Act3_78B1_15A_rEps3_1_rKdcKohc_5_Krl_365BL_Cdc0_0_25_Cohc0_0_25_80PT20000.mat');      
    end

    uBM=var.ModelOutputsStruct.structural.uBM;
    uTMR=var.ModelOutputsStruct.structural.uTMR;
    uTMB=var.ModelOutputsStruct.structural.uTMB;
    uDC=var.ModelOutputsStruct.structural.uDC;
    uStapes=var.ModelOutputsStruct.middleEar.uStapes;

    phiSV=var.ModelOutputsStruct.electrical.phiSV;
    phiSM=var.ModelOutputsStruct.electrical.phiSM;
    phiOHC=var.ModelOutputsStruct.electrical.phiOHC;
    phiST=var.ModelOutputsStruct.electrical.phiST;

    xMesh=var.ModelGeneratedVarsStruct.xMesh;
    T=var.ModelOutputsStruct.T;
    F0=var.StimulusOptionsStruct.toneBurst.F;
    [~,indexPeak]=max(max(abs(uBM),[],2));

    ModelOptionsStruct=var.ModelOptionsStruct;
    OptionalStruct=var.OptionalStruct;
    OptionalStruct.removeDCDOF=0; % not eliminated so that we call calculate DC disp.

    x=xMesh(indexPeak);
    fs=1/(T(2)-T(1));
   
    % smooth the data to visualize DC shifts   --------------------
    uBMsmoothed= smooth(uBM(indexPeak,:),1000);    
    uDCsmoothed= smooth(uDC(indexPeak,:),1000);    
    uTMBsmoothed= smooth(uTMB(indexPeak,:),1000);
    
    figure(100)
    subplot(1,Nmodel,iModel)
    plot(T*1e3,uStapes*1e7)
    xlabel('Time (ms)');ylabel('u_{stapes} (nm)')

    fig1=figure(1);
    subplot(3,Nmodel+1,2*(Nmodel+1)+iModel)
    plot(T*F0,uBM(indexPeak,:)*1e7,'k')
    hold on
    plot(T*F0,uBMsmoothed*1e7,'Color',blue_CB,'LineWidth',2)
    line([0 100],[0 0],'Color','k','LineStyle',':')
    xlabel('Time (cycles)')
    if iModel==1
        ylabel(['BM disp.',char(10),'(nm)'])
    end

    subplot(3,Nmodel+1,(Nmodel+1)+iModel)
    plot(T*F0,uDC(indexPeak,:)*1e7,'k')
    hold on
    plot(T*F0,uDCsmoothed*1e7,'m','LineWidth',2)
  
    line([0 100],[0 0],'Color','k','LineStyle',':')
    if iModel==1
        ylabel(['OHC/DC disp.',char(10),'(nm)']) 
    end

    subplot(3,Nmodel+1,iModel)
    plot(T*F0,uTMB(indexPeak,:)*1e7,'k')
    hold on
    plot(T*F0,uTMBsmoothed*1e7,'Color',green,'LineWidth',2)
    line([0 100],[0 0],'Color','k','LineStyle',':')
    if iModel==1
        ylabel(['TM disp.',char(10),'(nm)'])
    end
    
   
    % calculate and plot the FFTs of the displacements  ---------
    T0=1e-3;
    Tf=2.7e-3;
    [uBMfreq,freq]=computeFFT2_FlatTopWindow(uBM,T,T0,Tf);
    [uDCfreq,~]=computeFFT2_FlatTopWindow(uDC,T,T0,Tf);
    [uTMBfreq,~]=computeFFT2_FlatTopWindow(uTMB,T,T0,Tf);

    figure(2)
    subplot(3,(Nmodel+1),2*(Nmodel+1)+iModel)
    semilogy(freq/1000,abs(uBMfreq(indexPeak,:)*1e7))
    xlim([0 100])
    ylabel('BM disp. (nm)')
    xlabel('Frequency (kHz)')
      
    subplot(3,(Nmodel+1),iModel)
    semilogy(freq/1000,abs(uDCfreq(indexPeak,:)*1e7))
    xlim([0 100])
    ylabel('DC disp. (nm)')
    xlabel('Frequency (kHz)')

    subplot(3,(Nmodel+1),(Nmodel+1)+iModel)
    semilogy(freq/1000,abs(uTMBfreq(indexPeak,:)*1e7))
    xlim([0 100])
    ylabel('TM disp. (nm)')
    xlabel('Frequency (kHz)')


    % find the harmonics from the FFT
    index_f0=find(freq==F0);
    index_2f0=find(freq==2*F0);
    index_3f0=find(freq==3*F0);

    uTMBcomp(:,1)=uTMBfreq(:,1)/2;
    uTMBcomp(:,2)=uTMBfreq(:,index_f0);
    uTMBcomp(:,3)=uTMBfreq(:,index_2f0);
    uTMBcomp(:,4)=uTMBfreq(:,index_3f0);

    uDCcomp(:,1)=uDCfreq(:,1)/2;
    uDCcomp(:,2)=uDCfreq(:,index_f0);
    uDCcomp(:,3)=uDCfreq(:,index_2f0);
    uDCcomp(:,4)=uDCfreq(:,index_3f0);

    uBMcomp(:,1)=uBMfreq(:,1)/2;
    uBMcomp(:,2)=uBMfreq(:,index_f0);
    uBMcomp(:,3)=uBMfreq(:,index_2f0);
    uBMcomp(:,4)=uBMfreq(:,index_3f0);

    figure(3)
    subplot(3,(Nmodel+1),iModel)
    plot(xMesh,20*log10(abs(uTMBcomp(:,1)*1e7)))
    hold on
    plot(xMesh,20*log10(abs(uTMBcomp(:,2)*1e7)))
    plot(xMesh,20*log10(abs(uTMBcomp(:,3)*1e7)))
    plot(xMesh,20*log10(abs(uTMBcomp(:,4)*1e7)))
    xlabel('x (cm)');ylabel('TM disp. re 1nm')
    xlim([0 0.5])
    ylim([-70 30])
    legend('Tonic','F_0','2F_0','3F_0')

    subplot(3,(Nmodel+1),(Nmodel+1)+iModel)
    plot(xMesh,20*log10(abs(uDCcomp(:,1)*1e7)))
    hold on
    plot(xMesh,20*log10(abs(uDCcomp(:,2)*1e7)))
    plot(xMesh,20*log10(abs(uDCcomp(:,3)*1e7)))
    plot(xMesh,20*log10(abs(uDCcomp(:,4)*1e7)))
    xlabel('x (cm)');ylabel('DC disp. re 1nm')
    xlim([0 0.5])
    ylim([-70 30])

    subplot(3,(Nmodel+1),2*(Nmodel+1)+iModel)
    plot(xMesh,20*log10(abs(uBMcomp(:,1)*1e7)))
    hold on
    plot(xMesh,20*log10(abs(uBMcomp(:,2)*1e7)))
    plot(xMesh,20*log10(abs(uBMcomp(:,3)*1e7)))
    plot(xMesh,20*log10(abs(uBMcomp(:,4)*1e7)))
    xlabel('x (cm)');ylabel('BM disp. re 1nm')
    xlim([0 0.5])
    ylim([-70 30])


end


%% Plot Dewey 2021 waveforms
load('Dewey2021_data/WT_ex_Fig1D_E.mat')
f0= data.f1;
period = 1./f0;
if normalizedTime==1
    Texp=data.LIVE.TM.t./period;
else
    Texp=data.LIVE.TM.t*1e3;
end

figure(1)
subplot(3,(Nmodel+1),(Nmodel+1))
plot(Texp-nOffset,data.LIVE.RL.td,'k')
hold on
plot(Texp-nOffset,data.LIVE.RL.td_lp,'Color',green,'Linewidth',2)
line([0 100],[0 0],'Color','k','LineStyle',':')

subplot(3,(Nmodel+1),2*(Nmodel+1))
plot(Texp-nOffset,data.LIVE.OHC_DC.td,'k')
hold on
plot(Texp-nOffset,data.LIVE.OHC_DC.td_lp,'m','Linewidth',2)
line([0 100],[0 0],'Color','k','LineStyle',':')

subplot(3,(Nmodel+1),3*(Nmodel+1))
plot(Texp-nOffset,data.LIVE.BM.td,'k')
hold on
plot(Texp-nOffset,data.LIVE.BM.td_lp,'Color',blue_CB,'Linewidth',2)
line([0 100],[0 0],'Color','k','LineStyle',':')

xlabel('Time (ms)')

ylims = [-32 32];
ylimsModel=[-32 32];


for i = [(Nmodel+1) 2*(Nmodel+1) 3*(Nmodel+1)]
    subplot(3,(Nmodel+1),i)
   
    if i == 3
        title('Experiment','Fontsize',11)  
    end
end

for i = 1:3*(Nmodel+1)
    subplot(3,(Nmodel+1),i)
    set(gca,'Fontsize',FS)
    xlim([0 80])
     
    ylim(ylimsModel)
    yticks([-30:10:30])
    if or(or(i==1,i==4),i==7)
        yticklabels({'','-20','','0','','20',''});
    else
        yticklabels([]);
    end

    xticks([0 20 40 60 80 100])

    if i==1
        title('M1')
    elseif i==2
        title('M4')
    elseif i==(Nmodel+1)
        title('Experiments')
    end

    if or(i==(Nmodel+1)+1,or(i==(Nmodel+1)+2,i==(Nmodel+1)+3))
        text(4,-21,char('A'+i-1),'FontSize',14','FontWeight','bold')
    else
        text(4,21,char('A'+i-1),'FontSize',14','FontWeight','bold')
    end

    if i>=2*(Nmodel+1)+1
        subplot(3,(Nmodel+1),i)
        xlabel('Norm. Time (cycles)')
        xticklabels({'0','','40','','80',''})
        xtickangle(0)
    else
        xticklabels([]);
    end
   
end

set(fig1,'Units','inches','Position',[2 2 6.5 4])








